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1 Intoduction 



This document contains the documentation for the MATLAB toolbox EXPODE. The toolbox 
provides advanced exponential integration methods featuring five different integrator classes, 
evaluation of the Matrix functions directly and by a Krylov subspace method and an adaptive 
step size implementation for some of the integrators. EXPODE is as compatible as possible 
to MATLAB's internal ODE toolbox syntactically and semantically, such that a switch to this 
package is made easy. 

For the mathematical details on the implementation we refer to jl]. We first give a short 
quick start guide in the remainder of this section and continue with the full documentation 
in the next section. Note that some information might be duplicated. 

1.1 Installation and Requirements 

We will now describe the minimal requirements and the installation of the EXPODE toolbox. 
EXPODE runs on all recent and middle-aged computers. The performance highly depends on 
the problem and available hardware. The toolbox was tested on MATLAB versions down to 
MATLAB 7.2 (R2006a), released in 2006. Older versions might be albe to run it as well and we 
would appreciate feedback on success or failure. Versions prior to 7.0 cannot be compatible 
due to the lack of proper function handles. OCTAVE is also not capable of running EXPODE, 
since it currently does not support nested functions which are used regularly in the toolbox. 
At 

http: //www. am.uni-duesseldorf . de/en/Research/03_Sof tware .php 

you can download the packages containing the EXPODE toolbox. Two different versions are 
available, a package for users and an extended one for developers. The latter contains some 
additional tools helpful for extending EXPODE. Usually the user package should be sufficient. 

To install, just unpack the archive obtained from the above link. In a UNIX environment 
you can do this by typing 

$ tar xzvf expode-VERSION.tgz 

in the directory where you downloaded the package. This will give you an expode subdirec- 
tory. To make it available in MATLAB, just add the package's root to MATLAB's path and run 
the initPaths function with 

>> addpath /download/path/expode; 
>> initPaths; 

To make a permanent installation for the current user, put the above line into your startup .m 
file. See MATLAB's help for more information. 
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1.2 Quick Start 



To get a first impression of EXPODE we start witfi running some of the included examples. To 
access the examples we add the examples directory to the MATLAB path. Run the following 
lines 

» addpath /download/path/expode/examples; 
» [t, y] = HeatlD([] , [] , 'run'); 

to solve a Heat equation with a time-dependent source term in one dimension. Use 
» help HeatlD; 

to get some more information on the equation. The solution over time will be plotted in 
a mesh plot. All examples contained in the package can be run by simply calling them 
without arguments. Short information on the equations are contained in their helptexts. To 
experiment with solver it is convenient that we now run the example manually. 

» °/o paramerters 

» epsilon =0.1; gamma = 0.1; N = 100; 
» 

» °/o get initial conditions 

» [tspan, yO, options] = HeatlD([], [] , 'init', epsilon, gamma, N) ; 
» 

» °/o run the example 

>> [t, y] = expode (OHeatlD, tspan, yO, options); 

Now we can start playing with options and parameters. Switching to direct the solver for 
the matrix functions, we use 

» options = expset (options, 'MatrixFunctions ' , 'direct'); 

Other options are set similarly. To get interactive help for options, we use one of the 
integrator specific set commands, i.e. 

» options = exprbset (options, 'MatrixFunctions', 'direct'); 
» options = exprbset (options, 'MinStep', -1); 

This will do some checks on the values set for an option. The second line will throw an error, 
because step sizes must be positive. An overview of the available options for an integrator 
is available by calling the integrator info without arguments. More detailed information on 
a specific option can be shown with this command as well: 
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>> exprbinfo % prints all available options for exprb 

>> exprbinfo MinStep % prints helptext for MinStep option 

The following code will produce a logarithmic error vs. time step plot. We need to solve 
an equation where we can evaluate an exact solution. This is done by calling ode(t, [] , 
'exact'). We use the semil example here, see its helptext for information on the equation, 
'exprb' indicates to only plot results for the Rosenbrock-type methods, the " uses the 
direct solver for the matrix functions and N = 50 is the number of gridpoints for the spatial 
discretization of semil's underlying partial differential equation. 

» allMethodsCOsemil, 'exprb', ", [] , 50); 

The easiest way to implement an own differential equation is to start with some tutorial files 
in examples/Hello_World. Then you can modify example files in the examples directory. 
MinEx.m or Template. m contain a lot of comments for an easy start. MinEx.m is simpler, 
where as Template. m uses more advanced features. 



2 The expode Integrator 

The EXPODE integrator package solves a given ordinary differential equation of the form 

y' = F{t,y) = Ay + g{t,y) 

with a number of exponential time integration methods, where the second form of the equa- 
tion is only required for the semilinear solvers. For details on a specific method, see one of 
the later sections. The option |integrator| is used to determine the actual integrator to use, 



see sections E] and m for details on how to set und use it. 

After downloading and extracting the code package, you will get a directory called expode. 
This directory contains the expode integrator and the options helper expset. This is the ge- 
neral integrator, which calls the appropriate integrator as needed. Then for each supported 
integrator class there exists a dedicated integrator command (e.g. exprb), an options helper 
that understands the specific options for this integrator (e.g. exprbset) and an info com- 
mand (e.g. exprbinfo) to support the user. Furthermore, there is an analogon to MATLAB's 
deval, devalexp and the function initPaths to setup MATLAB's searchpath. Additionally 
there are some subdirectories. One of these is the example directory. The programs in that 
directory can be run directly and be used as a reference how to call the integrator. 

To use expode you either have to move the package's content to your MATLAB working direc- 
tory or use 

>> addpath /path/to/expode 

to tell MATLAB where to find the integrator. You can also use 
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» addpath /path/to/expode/examples 

to directly run one of the examples. 

expode's user interface is adapted to MATLAB's internal integrators. A call with all available 
arguments is of the following form: 

>> [t, y] = expode((8)ode, {@jac}, {tspan, yo, opts}, {varargin}) ; 

To keep the explanation well-arranged, the possible call combinations will be explained step 
by step. The simplest way to invoke the integrator is 

» [t, y] = expode((Q)ode, tspan, yo) ; 

where ode is a function to evaluate the differential equation (and its first derivative) at given 
data y and t. tspan = [to, T] defines the integration interval. You can also use tspan = 
[to , . . . , tfj] for some integrators instead. This will result in a solution evaluated at the 
given times instead of the ones chosen by the integrator internally. In that case, clearly we 
have t = tspan. This kind of solution is called dense output because one would typically use 
tspan = linspace(to, T, N + 1) with different output data than the step selection would 
produce. exp4 allows dense output directly as it has an own dense output formula. For all 
other integrators there is only a generic hermite interpolation formula which is not suitable 
for stiff problems. This is disabled by default, but can be switched on via the IDDGeneratof] 
option. Use this feature with care, yo is the initial condition to solve the equation with. 

The ode function has to be callable in the following way: 
» res = ode(t, y, {flag}); 

t and y are the time and phase space variables respectively, flag controls the output of the 
function. An omitted or empty flag should trigger the evaluation of the right hand side. 
If the integrator used is for semilinear problems (exprk and expmssemi) flag = 'linop' 
should trigger the evaluation of the linear part (the matrix). In case of the linearized integra- 
tors exprb, expms and exp4, flag can be given as flag = Jacobian' or flag = 'df _dt'. 
In the first case the derivative of the right hand side with respect to y - here called Jacobian, 
see the introduction - has to be evaluated at t and y. The last case is only needed if the diffe- 
rential equation is non-autonomous. Then ode should return the derivative of the right hand 
side with respect to the time variable t. If the differential equation is non-autonomous, you 
have to set the option [NonAutonomousI to 'on', see the options section. A typical function 
body for ode is presented for a linearized integrator below. 

function res = ode(t, y, flag) 

if nargin == 2 I I isempty (f lag) 
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flag = 

end 

switch flag 
case 

res = evaluation of the right hand side; 
case 'jacobian' 

res = evaluation of the Jacobian; 
case iinop 

res = evaluation of the linear part; 
case 'df_dt' 

res = evaluation of the derivative of the 
right hand side w.r.t. t; 

otherwise 

error ( 'Unknown flag: 7oS . ' , flag); 

end 



return res ; 

end 



If you want to use both integrator types, simply add cases for ' jacobian', 'df„dt' and 
'linop'. The otherwise case ist very helpful to find errors. All function handles used can 
be replaced by inlines or the functions' names as strings. All of these can be evaluated by 
the EXPODE integrators. 

It is possible to source out the evaluation of the Jacobian or the linear part into an additional 
function. In that case, expode has to be called this way: 



>> [t , y] = expode (@ode, @jac, tspan, yo) ; 



where function J = jac(t, y) is the new function for the Jacobian. The same applies to 
an evaluation function of the linear part, which will be called jac as well for simplicity. Note 
that for consistency reasons this function should allow parameters t and y as well, which 
can be ignored by the code. The handle to these functions can alternatively be given as an 
integrator option. See the I JacobianI option in section or the [LinOp] option in section \^?2\ 
for details. If these options are switched 'off' and I JacobianVI or |LinOpV| are 'on' instead. 



jac will be interpreted as a function to evaluate the product of the Jacobian or linear part 
with a vector. In that case, function res = jac(t, y, v) has to be callable. 

To use integrator options, invoke expode in the following way: 



>> [t, y] = expode (@ode, {@jac}, tspan, yo, opts); 

The construction of such an options object will be discussed in the Inext section! 
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If you want to hand over additional parameters to the ode function - stiffness parameters for 
springs or some dimensions for example - these parameters can be passed to the integrator 
as follows: 

>> [t, y] = expode(@ode, {@jac}, tspan, yo, opts, varargin) ; 

varargin can be an arbitrary number of additional argument separated by comma. These 
parameters will be passed to the ode function directly via 

» res = ode(t, y, flag, varargin{ : }) ; 

Please do not use a struct object as the first of the varargin arguments if you do not use 
integrator options. The integrator would confuse this argument with an options structure. 

If you still want to do so, you have to provide an empty vector ( [] ) as the options parameter. 

function res = ode(t, y, flag, varargin) 

if nargin == 2 I I isempty (f lag) I I "ischar (f lag) 
if "ischar (flag) 

varargin = {flag, varargin{ : }} ; 

end 

flag = ; 

end 

"L The rest can be the same as Ibef orel . 

end 

If varargin is used, it will also be passed to the jac function if it is available. It will be 
called with J = jac(t, y, varargin{:}) in that case. 

It is also possible to omit the integration interval, the initial condition, and the integrator 
options if the ode function is able to provide them. If varargin is used, you need to use an 
empty vector ( [] ) for the omitted parameters. Calling expode with 

» [t, y] = expode (@ode, {@jac}); 
or 

» [t, y] = expode (@ode, {@jac}, [], [] , [] , {varargin}); 
will invoke ode with the ' ini : flag: 

» [tspan, yo, opts] = ode([], [] , 'init', {varargin}); 
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The example code for the ode function has to be extended by another case for ' init ' , as 
can be seen here. 



case 'init' 

res{l} = [to, T] ; 

res-[2} = yo; 

opts = exprbset (optionl , valuel, option2, value2, 

°/„ see section [3], |exprbse"t 

res{3} = opts; 



To be able to evaluate the solution at an arbitrary time after the integration is finished, the 
call 

» sol = EXPODE(@ode, {@jac}, {tspan, yo, opts>, {varargin}) ; 



will generate a variable sol which can be used with the devalexp function via 



>> [ y, dydt ] = devalexp(sol , t) ; 

The argument t then is a vector of times in the integration interval. Then y contains the 
evaluation of the numerical solution and dydt its first derivative ot the times in t. Depending 
on the dense output generator used, sol can get quite large. For the exp4 integrator it has 
to save all inner stages for each integration step, that is a total eight vectors per step. Note 
that this only works if a dense output formula is available. 



3 Using Options 

For each of our integrators we have a command to set options. Exemplary we show their 
usage with exprbset, the one for the exprb integrator. It is an extension to MATLAB's 
odeset. It would have been preferable to use odeset directly, but it has a fixed set of 
available options and is not extensible to support a larger set. 

Both odeset and exprbset create a so called options structure or options object, which 
contains the options set by the user. These structures can be passed to the integrators. It 



has been explained in the previous section how to do that for expode 



exprbset is compatible with the option objects created by odeset, so one can create an 
options structure with odeset and then extend it with exprbset. Unfortunately, it is not 
possible the other way around, since odeset removes all options it does not know about. The 
first sequence is more important though. Usually, a program will be written using MATLAB's 
standard tools and then be extended to use other integrators. 

The exprbset function can be called analogously to odeset. The simplest way to use it is 
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» opts = exprbset (option , value); 



where option is an option's name and value the user's choice. The latter one can be scalar, 
vector valued, logical, a string, or a function, depending on the option chosen. All supported 
option types are listed later in this section. The info commands such as exprbinf o provided 
by the EXPODE package can print options supported by the EXPODE integrators. 

It is also possible to pass more option- value-pairs at once and extend already existing option 
objects: 



» opts = exprbset (optionl , valuel, option2, value2, ...); 
» opts = exprbset (opts , optionS, valueS, ...); 



For a more comfortable usage, you can use the integrator specific set commands (exprbset, 
exprkset, etc.) to check if the provided options are valid, expset itself cannot do this, since 
it does not know which integrator will be used. You can tell it by using 



» opts = exprbset (intname , {opts}, optionl, valuel, option2, value2, ...); 



where intname is the name of the integrator, e.g. 'exprb', 'exprk'. 

In the remaining part of this section, the supported value types for options will be discussed. 
All of the listed types can be checked for validity by the set commands. Which type is 
allowed for which option is described in the Inext sections! The info commands in MATLAB 
can provide this information. 



» exprbinf (-[optname}-) ; 



will display a list of all options supported by exprb if optname is omitted. If it is given, 
more detailed help on the specific option will be printed. Using optname = ' - ' will give 
this detailed help on all options. Warning: the output is quite long, it should be used with 



» more on; exprbinf o( ); more off; 



The output of the info command contains the option's type in the first line in squared 
brackets. For the lAbsToll option this line reads 

AbsTol - Absolute error tolerance [ positive scalar I positive vector {le-06}- ]. 
The vertical bar separates alternative types, so the type is |positive|lscalarl or [positive" 



Ivectori The default value is set in braces after the type. The two types are each a combi- 
nation of two other types: |positiv"e] and Iscalarl or Ivectori respectively. Which types are 
compatible will be explained below. Another example for types is given by the IJacobianVI 
option. Here the exprbinf o command prints 
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JacobianV 



[ function.handle I {'off'} I 'on' ]. 



The allowed types in this case are IbooleanI or [f unction_handle| For types Ibooleanl and its 
generalization llistl a list of accepted values will be stated. For a Ibooleanl this list is fixed 
and consists of the values (true) and 'off (off), for the llistl type the list depends 

on the option. The listed values are always quoted strings. Now all option types will be 
discussed in detail. 



Type scalar. This type can be any numerical type that is supported by MATLAB. scalar 
can be combined with jinteger"} jpositivej |non-negative[ jnegative] and jnon-positive 



Type vector. This type can be a vector of any numerical type that is supported by 
MATLAB. Iscalarl is the special case of this type with length one, so any scalar is a vector 
as well. Vector can be combined with 



integer positive, non-negative, negative 



and 



non-positive 



Type matrix. This type can be a matrix of any numerical type that is supported by 
MATLAB. Iscalarl and Ivectorl are special cases of this type, so any scalar and vector is a 
matrix as well. Matrix can be combined with|integer[ [pbsitive[ |non-negative[ |negativ"e 



and jnon-positive 



Type integer. This type can be any integer number. This does not enforce the use 
of MATLAB's intS, intl6 or similar types, value is accepted as an integer as long as 
the expression (value - round (value) ) == ist true, integer can be combined with 
Ivector^ Iscalarl jpositivej jnon-negativej, jnegative] and jnon-positive 



Type positive. This type can be any numerical type that is supported by MATLAB (scalar, 
vector or matrix) where all components are strictly positive, positive can be combined with 
Imatrixj Ivectorl Iscalarl and jinteger" 



Type non-negative. This type allows any numerical MATLAB type (scalar, vector or 
matrix) where all components are not negative (positive or zero), non-negative can be 
combined with Imatrixj Ivectorl Iscalarl and 



mtegerj 



Type negative. This type can be any numerical type that is supported by MATLAB (scalar, 
vector or matrix) where all components are strictly negative, negative can be combined 
with Imatrixj Ivectorl Iscalarl and 



integer 



Type non-positive. This type allows any numerical MATLAB type (scalar, vector or 
matrix) where all components are not positive (negative or zero), non-positive can be 
combined with Imatrixj Ivectorl Iscalarl and jinteger 
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Type index. This type is an abbreviation for |positive|[integer|[scalar[ index can 
therefore not be combined with any other type. 



Type indices. This type is an abbreviation for |positive|[integer|lvectorl indices 
therefore cannot be combined with any other type. 



Type boolean. The values of this type can be given in different forms. Either use 'on' 
and 'off' or use the numerical values true and false, boolean is not combinable with 
any other type. Alternative string values are 'yes' and 'no' or 'true- and false'. 



Type list. The values for this type can be given in different forms. Either use one of the 
available values printed by the info command or use the number of the list entry, start coun- 
ting from 0. Example: The IMatrixFunctionsI option's info states (ignore f unction_handle 
here, it will be discussed below) 

MatrixFunctions - ... [ {'direct'} I 'arnoldi' I f unction_handle ] 

which permits to set the values 'direct' and 'arnoldi', in that order. Therefore 'direct' 
has number and 'arnoldi' has number 1. Some special list type options define spe- 
cific numeric values, which will be stated in parentheses behind the string value. See, for 
instance, option [Order] of the exprb integrator. It is possible to set the value for the option 
IMatrixFunctionsI to use Arnoldi's method in the following ways: 



>> opts = exprbset ( 'MatrixFunctions ' , 'arnoldi'); % or equivalently 
>> opts = exprbset ( 'MatrixFunctions ' , 1); 

list cannot be combined with any other type. 



Type text. This type can be any MATLAB char array (string), text cannot be combined 
with any other type. 



Type struct. This type can be any MATLAB struct object, struct cannot be combined 
with any other type. 



Type f unction_hajidle. This type can be either a MATLAB function handle, a MATLAB inline 
function or a char array (string). A function handle to a function f unc is generated by Of unc. 
An inline function is generated by MATLAB's inline command. If the argument is a string, 
it has to represent a function's name: use 'func for the function func. function_handle 
cannot be combined with any other type. 
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4 Available Options 



In the following sections, the options available for the EXPODE integrators will be discussed. 
The first option discussed will be an exception to this, as it is used to choose the actual 
integrator when calling expode directly. The rest of this section will be split into several 
subsections, starting with options common to all integrators, followed by options for semi- 
linear, then linearized integrators. The last two subsections are dedicated to constant and 
variable step size integrators. Let's start with the integrator option. 



Option Integrator. Integrator to use 

Permitted values: {'exprb'}, 'expms', ' expmssemi ' , 'exprk' or 'exp4' 

Select the integrator to do the actual work. Simply type the integrator's name in the MATLAB 
prompt to get more information. Instead of using this options it is preferred to run the 
appropriate integrator directly. This option mainly exists for automation purposes in scripts. 



4.1 Options Common to All Integrators 

Here we describe options common to all integrators. The options will be sorted into several 
groups. First some options for controlling the integration process, then some properties for 
the ODE followed by output control and debugging options. 



Option AbsTol. Absolute error tolerance 

Permitted type: |positive]lscalarl or |positive|lvectorl -rie-06> 

A scalar tolerance applies to all components of the solution vector. Elements of a vector 
of tolerances apply to corresponding components of the solution vector. lAbsToll defaults to 
10~^ in all solvers. This also controls the stopping criterion in the Krylov process if matrix 
functions are evaluated this way. 

See also: IRelToll and INormControI] 



Option RelTol. Relative error tolerance 
Permitted type: |positive]lscalarl lO . OOll- 

This scalar applies to all components of the solution vector, and it defaults to 10"'^ (0.1% 







accuracy) in all solvers. The estimated error in each integration step satisfies \\err\\ < 1 in a 
norm scaled with RelTol •max(abs(yn(i), abs(yn_i(z))) -l-AbsTol(z) in each component, where 
is the numerical solution at the current, yn-i the one at the previous time step. This also 
controls the stopping criterion in the Krylov process if matrix functions are evaluated this 
way. 

See also: lAbsToll and INormControI] 



13 



Option NormControl. Control only 2-norm of the error 
Permitted values: {'off'} or 'or 

Set tliis property to ' on ' to request tliat tlie solver controls the error in each integration 
step to meet ||err|| < (1/n) ■ max{RelTol ■ \\y\\, AbsTol). By default the solvers use a more 
stringent componentwise error control. 

See also: lAbsToll and IRelToll 



Option MatrixFunctions. Evaluation method for the matrix functions 
Permitted values: {'direct'}, 'arnoldi or |f unction_handle| 

Set the method to evaluate the product of the matrix function with vectors. There are 
two built-in methods. The default direct ' method uses diagonalization of the Jacobian 
to do this. If the Jacobian is too large, this will be too expensive computionally - use 
'arnoldi ' in this case. All matrix functions will then be approximated in Krylov subspaces, 
the Arnoldi method is then used to compute a nested orthonormal basis of that space. 
Another alternative is to provide a custom function to compute these results. See the section 
[To] how this method has to work. 



See also: IJacobiant IJacobianVL |LinOp[ |LinOpV| and [KrylovTest Index 



Option KrylovTestlndex. Dimensions of the Krylov subspaces to test the residual 
Permitted type: lindexl or fvectorllofllindiciesl { [ 1 2 3 4 6 8 11 15 20 27 36 46 57 
70 85 100 ]} 

Sequence of ascending integers which indicate the dimensions of the Krylov subspaces where 
the residual is tested. If the KrylovMaxDim option is set to zero, the largest number is 
the maximum dimension, otherwise dimensions greater than specified there will be ignored. 
For KrylovMaxDim itself, the residual will also be tested. If the Krylov process does not 
terminate within this range, the time step size is reduced. This option only applies if 
MatrixFunctions is set to ' arnoldi ' . 



See also: IMatrixFunctionsl |KrylovMaxDini| and [KrylovAbort 



Properties of the differential equation can be set by the following options. 



Option NonAutonomous. Specifies whether the ODE is autonomous or not 
Permitted values: {'off'} or 'on' 

Set to „ if the differential equation is not autonomous. This option is only required for 

integrators using the Jacobian. It prints a message at integrator startup for all integrators 
though. 



Option Complex. Solution is complex 
Permitted values: or { } 
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Set to 'off' if the solution does not have complex components. This will set all possibly 
numerically generated imaginary parts of the solution to zero. 



Option Structure. Structure of the Jacobian/linear part 

Permitted values: {'none'}, 'normal', 'symmetric', ' skewsymmetric ' or 'diagonal' 

If the Jacobian (or linear part) has a special structure, this can be exploited in the matrix 
function evaluators. Available properties are none - if there is no special structure, diagonal 

- if the matrix is diagonal, symmetric - for symmetric or Hermitian matrices, skewsymmetric 

- for skew-symmetric or skew-Hermitian matrices and normal - for normal matrices, that 
are none of the three previous. 



Aliases: |Symmetr"y 



Option GFcn. Evaluation of the nonlinear part of the ODE, for semilinear problems 
Permitted type and values: |f unction_handlel 'off or {'od } 

This option will only be used, if ISemiliiil is ' on ' . Set this to ' on ' if the ODE file can 
evaluate the nonlinear part when used with flag 'gfun'. Set to a function handle if an 
external function can evaluate it. Setting to will result in reverting to the standard 

solving method, ignoring the special structure of the ODE. 

See also: ISemiliiil 

Now we will discuss how to control the output of the integrator. 
Option DOGenerator. Dense Output generator 

Permitted type and values: |f unction_handlel { default '} or 'hermite' 

Set the method to generate dense output, 'default' is the integrator's default method, 
which is none for the most integrators, 'hermite' uses hermite interpolation and is usa- 
ble with all integrators, but is not applicable for stiff problems! Use it for testing only. 
Integrators having a designated dense output formula will override and extend this option. 

To obtain solutions at specific times tO,tl,...,tf inal (increasing or decreasing) use tspan = 
[to tl ... tf inal] when calling 

>> [tout, yout]=expode (ode, tspan, yO, . . . ) ; 

or 

>> [tout ,yout] =expode (ode ,jac , tspan, yO, . . . ) ; 
>> sol=expode(ode, . . . ) ; 

will generate a variable sol which can be used with the devalexp function via 
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» [y ,dy] =devalexp(sol ,t) ; 

with a vector of times t. You will get the numerical solution and its first derivative evaluated 
at the times in t. 



Option Refine. Output refinement factor 
Permitted type: |positive||integer]lscaliLr] -Til- 



This property increases the number of output points by the specified factor producing smoo- 
ther output. Refine defaults to 1. Refine does not apply if length(tspan) > 2. To use this 
feature, you need a dense output generator. 

See also: IDOGeneratorl 



Option OutputFcn. Output function called after each time step 
Permitted type and value: |f unction_handle] or { ' of f ' } 

If a function handle is given, it is called after each time step. The output function has to 
understand the following calls: 



» outputFunctionC [to ,tf inal] ,y, 'init' ,varargin) 
» outputFunction(t ,y , ,varargin) 



The first call should initialize the output function. The second call is executed at every time 
step - either the ones selected by the solver's error estimator, at the refined steps or at the 
steps definded by tspan argument at the solver's call, to, tf inal and t are scalar times, 
to and tf inal limit the integration interval and t is the time of the current step, y is the 
solution at to on the 'init call and the solution at t on init=". varargin is an argument 
which will be forwarded from the solver call. 



See also: [OutputSeT 



Option OutputSel. Indices of the solution given to the output function 
Permitted type: Ivectorllof IlindicesI { [] } 

Only used if an output function is used. If this argument is an empty vector, all indices of the 
solution are passed to the output function. Otherwise, y(OutputSel) is passed. OutputSel = 
1 : length(y) will give the same result as OutputSel = []. 

See also: |OutputF"cn 

The last set of options in this section deals with debugging and logging. 



Option Stats. Display status messages 

Permitted values: 'silent', 'off', {'on } or 'verbose' 
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Set the level of status messages printed by the integrator. If IJacobianStat'sl [StepStats 

' auto 



and |MatrixFunctionLog| are set to 'auto , they will be activated only when IStatsI is set 
to 'verbose' and be deactivated on all other settings. Setting this option to 'silent' will 
also disable warnings. 

See also: IJacobianStatsL |StepStats| and IMatrixFunctionStatsI 



Option MatrixFunctionStats. Display status messages for matrix function evaluations 
Permitted values: oti ', -on' or {'auto'} 

Control, whether or not to display status messages for the matrix function evaluations. When 
set to auto, this will be activated when [St at si is set to 'verbose'. 

See also: IStatsI 



Option Waitbar. Show a waitbar to display progress 
Permitted values: { }, ou , 'text' or 'both' 

Set to on' to display the current progress after each time step. Set to 'text' if you want 
to display the progess on the MATLAB prompt. Set to 'both' if you want both kind of 
displays. 

See also: IStatsI 



Option ClearlnternalData. Clear internal integrator data 
Permitted values: ' of f ' or { ' on ' } 

Set to 'of f ' if you want to keep the integrator's internal data. This data is available in the 
global eD variable. 



4.2 Options for Semilinear Integrators 

In this subsection we will discuss options common to the semilinear integrators exprk and 
expmssemi. To evaluate the linear operator we have the following two options. 



Option LinOp. Evaluation function for the ODE's right hand side's linear part 
Permitted type and values: |f unction_handlel ' of f ' , { ' on ' } or Imatrixl 

The integrator needs to compute matrix functions with the linear part of the right hand 
side of the ODE. Therefore, either the linear part has to be evaluated, or the evaluation 
of the linear part times vector has to be provided - this only works when using Krylov 
approximations to the matrix functions. If this option is set to on' (default), the ode 
function will be called with 



>> lin=ode(t ,y, 'linop' ,varargin) ; 
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If it is a function handle, the call 



» lp=opts .LinOp; lin=lp(t ,y ,varargin) ; 

will be executed. Note that the linear part will only be evaluated once. 
See also: |LinOpV| and IMatrixFunctionsI 



Aliases: IJacobianl 

Option LinOpV. Evaluation function for the ODE's right hand side's linear part times vector 
Permitted type and values: |f unction_handlel {'off'} or 'on' 

When using Krylov approximations to the matrix functions of the linear part of the right 
hand side times vectors, this option can be used to provide the result of the computation of 
the linear part multiplied by a vector. Sometimes it is faster to compute this directly instead 
of computing the full linear part. If this option is set to 'on' (default), the ODE function 
will be called with 

» res=ode (t , y , 1 inpop_v ' , v , varargin) ; 

If it is a function handle, the call 

» lpv=opts .LinOpV; res=lpv(t ,y,v, varargin) ; 
will be executed. 

See also: |LinOp| and IMatrixFunctionsI 



Aliases: IJacobianVI 

Additionally you can control the logging output with this option. 

Option LinOpStats. Display status messages for the evaluation of the linear part 
Permitted values: off, 'on' or {'auto'} 

Control, whether or not to display status messages for the evaluation of the linear part. 
When set to auto, this will be activated when [Stat si is set to 'verbose'. If enabled, this 
options prints 

» Evaluating linear part . . . done 

when the linear part is evaluated. If option [Stat si set to 'verbose' the Eigenvalues of the 
linear part will be plotted additionally. 

See also: IStatsI 

Aliases: IJacobianStatsI 
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4.3 Options for Linearized Integrators 

In this subsection we will discuss options common to the linearized integrators exprb, expms 
and exp4. Again the options are split into smaller groups to keep the information clearly pre- 
sented. The first couple of options deal with the evaluation of the needed data. Additionally 
we have one additional option for the equation's properties and one for logging. 



Option Jacobian. Evaluation function for the ODE's right hand side's Jacobian 
Permitted type and values: |f unction_handlel ' of f ' , { ' on ' } or Imatrixl 

The integrator needs to compute matrix functions with the Jacobian of the right hand 
side of the ODE. Therefore, either the Jacobian has to be evaluated, or the evaluation of the 
Jacobian times vector has to be provided - this only works when using Krylov approximations 
to the matrix functions. If this option is set to (default), the ode function will be called 
with 

» j=ode(t ,y , 'jacobian' ,varargin) ; 
If it is a function handle, the call 

>> jac=opts . Jacobian; j=jac(t,y,varargin) ; 

will be executed. If the Jacobian is constant, set this option to the constant matrix. 
See also: IJacobianVI and IMatrixFunctions] 



Aliases: [LinOp 



Option JacobianV. Evaluation function for the ODE's right hand side's Jacobian times 
vector 

Permitted type and values: |f unction_handlel { ' } or 'on' 

When using Krylov approximations to the matrix functions of the Jacobian of the right hand 
side times vectors, this option can be used to provide the result of the computation of the 
Jacobian multiplied by a vector. Sometimes it is faster to compute this directly instead of 
computing the full Jacobian. If this option is set to 'on' (default), the ODE function will 
be called with 



>> res=ode(t ,y, ' jacobian_v , v , varargin) ; 



If it is a function handle, the call 



» jacv=opts . JacobianV; res= jacv(t,y,v, varargin) ; 
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will be executed. 

See also: IJacobianl and iMatrixFunctions] 
Aliases: [LinOpV 



Option GJacobian. Evaluation of the Jacobian of the nonlinear part of the ODE, for 
semilinear problems 

Permitted type and values: |f unction_handlel 'off or {'on'} 

This option will only be used if ISemilTnl is 'on'. Set this to 'on' if the ODE file can 
evaluate the Jacobian of the nonlinear part when used with flag 'dg_dy'. Set to a function 
handle if an external function can evaluate it. Setting this to you will need to use 

Krylov approximations to the matrix functions and set the IGJacobianVI to something else 
than 'off. 

See also: ISemiliiil IGJacobieinVI and IMatrixFunctions] 



Option GJacobianV. Evaluation of the Jacobian of the nonlinear part of the ODE times 
vector, for semilinear problems 

Permitted type and values: |f unction_handlel { } or 'on' 

This option will only be used if ISemilinl is 'on' and MatrixFunctions is not direct'. 
Set this to if the ODE file can evaluate the product of the Jacobian of the nonlinear 

part with a vector when used with flag 'dg_dy_v'. Set to a function handle if an external 
function can evaluate it. Setting to ~;ff will use IGJacobianl 

See also: ISemiliiit IGJacobianl and [MatrixFunctions] 



The differential equation can now have this additional property. 



Option Semilin. Specifies whether the ODE is semilinear 
Permitted values: {'off } or 'on' 

Set to , if the differential equation has the form 

y' = Ay + g{t,y) 

Then add the flags 'gfun' and dg_dy' to your ODE file to evaluate the non-linear part 
g and its derivative with respect to y. Alternatively, you can use the IGFunI and IGJacobianl 
options to set these. If you use Krylov approximations to the matrix exponentials, you 
can also add the flag 'dp; dv v' to evaluate the Jacobian of g times a vector or use the 
IGJacobianVI option to do so. 

See also: IGFcnl IGJacobianl IGJacobianVI IMatrixFunctions] and IJacobianVI 



Please note that for the semilinear integrators in the next section this option is always 
assumed ' on ' . 
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The last option extends the logging capabilities. 



Option JacobianStats. Display status messages for the evaluation of the Jacobian 
Permitted values: off, 'on' or {'auto'} 

Control, whether or not to display status messages for the evaluation of the Jacobian. When 
set to auto, this will be activated when IStatsI is set to 'verbose'. If enabled, this options 
prints 

>> Evaluating Jacobian . . . done 

each time the Jacobian is evaluated. If option [StatsI set to 'verbose' the Eigenvalues of 
the Jacobian will be plotted additionally. 

See also: IStatsI 



4.4 Options for Constant Step Size Integrators 

Now we can move on to integrators without step size control, that currently are exprk, expms 
and expms semi. Since there is not so much to control, we only have one option. 



Option StepSize. Step size to use 
Permitted type: |non-negative|[scalarl {.0} 

Stepsize to use. If set to (default), the integrator will choose a default step size of (tf inal — 
to)/100. 

Aliases: [InitialStep 



4.5 Options for Variable Step Size Integrators 

The last general set of options is for the variable step size integrators exprb and exp4. Here 
we have the following options. 



Option hConstant. Use constant step size 
Permitted values: { ' of f ' } or ' on ' 

Uses constant step size in the integrator. In case IhConstantI is set to 'on', consider setting 
the |lnitialStep| option to specify the stepsize. 



See also: [InitialStep 



Option InitialStep. Initial step size to use 
Permitted type: |non-negative|[scalarl {.0} 
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Initial step size to use. If set to 0, the integrator will choose a default step size of (tf inal — 
to)/100. If IhConstantI is 'on' then the given step size here will be used for the entire 
integration process. 

See also: IhConstantI [MinStep] and [MaxStep 



Aliases: [StepSize 



Option MaxStep. Maximal step size to use 
Permitted type: |non-negative|[scalarl {.0} 



Set the maximal step size to use in the integration process here. Setting [MaxStep] to 
(default) will result in using (tfinal — to)/10, so using at least ten integration steps until 
the integrator is finished. This option will be ignored when using constant step size. 



See also: |MinStep"| |InitialStep| and IhConstantI 



Option MinStep. Minimal step size to use 
Permitted type: |non-negative|[scalarl {.0} 

Set the minimal step size to use in the integration process here. Setting [MinStep] to 
(default) will result in using eps(t) at integration time t, so that t + h is at least different 
from t. This option will be ignored when using constant step size. 

See also: [MaxStep] [InitialStep[ and IhConstantI 



And we have another option for logging again. 



Option StepStats. Display status messages for step size related events 
Permitted values: 'off, 'on' or {'auto } 

Control, whether or not to display step size related status messages. When set to auto, this 
will be activated when [St at si is set to 'verbose'. You will be informed about step size 
reductions and step rejections if this option is activated. 

See also: IStatsI 



5 exprk — The Exponential Runge-Kutta Integrator 

Since we are now familiar with the common options, it is now time to move to the actual 
integrators. This section will start with the exponential Runge-Kutta methods. EXPODE 
only supports explicit schemes. A good overview over this class can be found in [3]. exprk 
is a semilinear solver and only supports constant step size integrators. This can easily be 
extended to variable step size if appropriate error estimators are available for the schemes. 

To call exprk, use the following syntax: 
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>> [t, y] = exprk(@ode, {tspan, yo, opts}-, {varargin}) ; 
See section |2] for details on the arguments. 

The exprk integrator supports a number of Runge-Kutta schemes. See [3] for a reference. 
The scheme to use can be selected via the ISchemel option, see below. It is possible to eit- 
her select one of the predefined schemes or even provide your own. To do so, you have to 
create a scheme with the rkScheme function from the exprk subdirectory of the EXPODE 
distribution. rkScheme has quite a detailed help built-in (help rkScheme in MATLAB). Use 
initPaths ( ' exprk ' ) to put the exprk directory into your MATLAB-Path. To better un- 
derstand the usage of this tool, its usage is shown by constructing Krodstad's [7] scheme 
here. 



% Make sure tp have /path/to/expode in your MATLAB path 

initPaths ( ' exprk ' ) ; 

°/o Create an empty 4-node scheme 

sc = rkScheme (4, 4); 

% Set the nodes 

sec = [ 1/2 1/2 1 ] ; 

7o U2 = u„ + 1/2 ipi{-C2 h A) (g(Ui) - A u„)) 
sc = rkScheme (sc, 2, 1, 1/2); 

7o Us = u„ + (1/2 (^i(-C3 h A) - v?2(-C3 h A)) (g(Ui) - A u„) + ... 
I V52(-C3 h A) (g(U2) - A u„) 

sc = rkScheme (sc, 3, 1, [ 1/2, -1 ]); 
sc = rkScheme (sc, 3, 2, [0, 1 ]); 
% Stage 4 

sc = rkScheme (sc, 4, 1, [1, -2 ] ) ; 
sc = rkScheme (sc, 4, 3, [0, 2 ]); 
% Outer stage 

sc = rkScheme (sc, 1, 'b', [1, -3, 4 ]); 
sc = rkScheme (sc, 2, 'b', [ 0, 2, -4 ] ) ; 
sc = rkScheme (sc, 3, 'b', [0, 2, -4 ] ) ; 
sc = rkScheme(sc, 4, 'b', [0, -1, 4 ]); 



For reference on more complex schemes that use (pk{ci) for stage Uij or other functions than 
the yp's, see the function getRKSchemeFromOption found in the exprk directory as well. 

In addition to the icommon^ constant step size and lsemilinearl options. there are two additional 
options to determine the Runge-Kutte scheme to use, which will be explained now. 



Option Scheme. Select Runge-Kutte scheme 

Permitted values: 'Euler', 'StrehmelWeinerA' , 'StrehmelWeinerB', 'HeunA', 'HeunB', 
'CoxMatthews', {'Krogstad'}, 'HochbruckOstermann' or [struct] 
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Select the exponential Runge-Kutta coefficients here. The exponential Euler has no internal 
stage and is ffist-order convergent. StrehmelWeinerA and StrehmelWeinerB are of second 
order, where StrehmelWeinerA is generally preferable due to some order-loss on certain 
properties of the nonlinear part of the differential equation. HeunA and HeunB are both of 
order three. CoxMatthews, Krogstad are schemes of classical order four, where HochbruckO- 
stermann has full order four. You can also supply you own scheme, generate one with the 
rkscheme helper function from the exprk subdirectory. 

See also: IStrehmelWeinerParameterl 



Option Parameters. Select the free parameter (s) for some of the Runge-Kutta schemes 
Permitted type: Iscalarl or Ivectdrl 

Select the free parameter in the StrehmelWeinerA, StrehmelWeinerB, HeunA and HeunB 
schemes. In particular, the ffist parameter is always the second node for the scheme. In case 
of HeunB, the second parameter is the additional free parameter gamma. 



6 exprb — The Exponential Rosenbrock-type Integrator 

Now we turn our attention to exponential Rosenbrock-type methods [B]. They belong to the 
class of linearized integrators and have an embedded error estimator, so have an adaptive 
time stepping. 

To call exprb, use the following syntax: 



» [t, y] = exprb (@ode, {@jac}, {tspan, yo, opts}, {varargin}) ; 



See section [2] for details on the arguments. 

Due to their simpler order conditions there are much less schemes. In addition to the 
icommonl [variable step size] and IlinearizationI options there are two additional options. 



Option Order. Integrator order to use 

Permitted values: (2), 'three' (3) or {'four'} (4) 

Order to use for the integrator. When setting to two, the exponential Rosenbrock-Euler 
method will be used. The exprb32 method is the order three method with two inner stages 
using the second order exponential Rosenbrock-Euler method as error estimator. For order 
four, the exprb43 method will be used. This one has three inner stages and uses a three 
staged order three error estimator. 



Option ErrorEstimate. Error Estimator Parameters 
Permitted type: [vector] { [ -2 ]} 
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The order four integrator uses an error estimator that can be tweaked with two parameters 
a and b. These two parameters can be specified here with lErrorEstimatel = [a, b] . Choose 
b 7^ 1 — 6a, otherwise the error estimator gets close to weak order four instead of three. 



7 expmssemi — The Exponential Multistep Integrator 

This section discusses the exponential multistep integrator expmssemi, cf. and [1]. As 
with exponential Runge-Kutta methods it's a semilinear solver and only supports constant 
step sizes. In contrast to these methods, it is much harder to generalize them to variable 
step sizes due to their construction. 

To call expmssemi, use the following syntax: 



» [t, y] = expmssemi (@ode , {tspan, yo, opts}, {varargin}) ; 



See section |2] for details on the arguments. 

In addition to the icommon^ constant step size and Isemilinearl options there is one additional 
option to determine the step count. 



Option kStep. Number of old time steps to use 

Permitted values: (1), (2), 'three' (3), {'four'} (4), 'five (5) or 'six' 

(6) 

Number of old time steps the method should use. One gives the exponential Euler method, 
which is first order convergent. The number of steps is automatically the order of the scheme 
as well. 



Option StartupSteps. Source of the startup steps needed for the multi step method 
Permitted values: {'Fixpoint'}, 'Exact' or 'ExpRK' 

Select the computation method for the startup steps for the multistep method. 'Fixpoint' 
uses a fixedpoint iteration to compute the steps. 'Exact' queries the ODE with the flag 
'exact'. 'ExpRK' uses an exponential Runge-Kutta scheme of the appropriate order. Note 
that only methods up to order five are possible this way. 



8 expms — The Exponential Linearized Multistep Integra- 
tor 

This section discusses the exponential linearized multistep integrator expms. It is a linearized 
variant of the expmssemi and only supports constant step sizes. We have implemented the 
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general scheme introduced by Hochbruck and Ostermann in [S] and a scheme proposed by 
Tokman in [9]. 

To call expms, use the following syntax: 



>> [t, y] = expms (@ode , {@jac}, {tspan, yo, opts}, {varargin}) ; 



See section [2] for details on the arguments. 

In addition to the icommon^ constant step size and llinearizat ioni opt ions there is one additional 
option to determine the integration scheme. 



Option kStep. Number of old time steps to use 

Permitted values: (1), (2), three' (3), {'four } (4), 'five' (5) or 'Tokman' 

(6) 

Number of old time steps the method should use. One gives the exponential Rosenbrock- 
Euler method, which is second order convergent. The number of steps plus one is auto- 
matically the order of the scheme as well. The 'Tokman' two step scheme was proposed 
earlier than the more general construction by Hochbruck and Ostermann and is third order 
convergent. 



Option StartupSteps. Source of the startup steps needed for the multi step method 
Permitted values: { iixpoint'}, 'Exact' or 'ExpRB' 

Select the computation method for the startup steps for the multistep method. 'Fixpoint ' 
uses a fixedpoint iteration to compute the steps, 'p^-act' queries the ODE with the flag 
'exact'. 'ExpRB' uses an exponential Rosenbrock-type scheme of the appropriate order. 
Note that only methods up to order five are possible this way. 



9 exp4 - Exponential Integrator of Order Four 

This section is dedicated to the exp4 integrator [2]. It belongs to the class of linearized 
integrators and has adaptive time stepping. 

To call exp4, use the following syntax: 

» [t, y] = exp4(@ode, {@jac}, {tspan, yo, opts}, {varargin}) ; 

See section [2] for details on the arguments. 

Since exp4 has its own dense output formula, it overrides the IDOGeneratorl option. 
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Option DOGenerator. Dense Output generator 

Permitted type and values: |f unction_handlel { exp4'}, 'hermite' or 'none' 

Set tlie metliod to generate dense output. Default is 'exp4', exp4's own dense output 
formula. It can be overriden to use hermite interpolation by setting this option to 'hermite ' . 
Hermite interpolation is not applicable for stiff problems. Use it for testing only. Additionally 
the dense output generator can be switched off. 

10 Matrix Functions 

This section deals with custom evaluation functions for the product of matrix functions with 
vectors. 

expode has two built-in methods to compute these evaluations: directly by diagonalisation 
and using a Krylov subspace method. In some situations, it can be useful to provide a 
custom function that does this job. If the Jacobian has a special structure that makes it 
possible to compute matrix exponentials differently even in large dimensions, this would be 
such a case. For this reason, it is possible to hook an arbitrary function into the integration 
process. See the IMatrixFunctionsI option on how to accomplish that. 

The two internal functions are implemented the same way as an external method would have 
to work. They can be used for reference and can be found in the EXPODE distribution at 
matPun/matPunDirect .m and matPun/matPunKrylov .m. The former one is much easier to 
understand due to its size and complexity. 

It will be neccessary to use the global eD, jac and gjac variables which contain informa- 
tion shared between the various EXPODE functions. See subsection 111.41 for details on these 
variables. The parts important for this context will be explained when they appear. 

All logging output of the function should be printed with the eD . log .matPunLog function. 
This way it can easily be redirected or disabled by the integrator. 

The custom function will be called matPun in the remainder of this section. It will be 
invoked in different ways by the EXPODE integrators. The function's signature needs to have 
the following form: 

function [ h, varargout ] = . . . 

raatPunCjob, t, y, h, flag, v, reusable, reuse, facs) 
global eD; 
= eD . int . o ; 
funs = eD . functions ; 

The flag variable has two different roles. It can either request a specific action to be 
performed or it can inform the method about the meaning of the currently evaluated matrix 
function. The request-type flags will be discussed first. 
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»matFun([], [] , [] , [] , 'init'); 



should initialize the function globally. It has to state whether or not it requires to explicitly 
evaluate the Jacobian of the right hand side and - if used - the one of the non-linear part g 
for semilinear problems (see option ISemilin|) . This is done by setting 

eD .matFun .needJacExplicit = true; 
eD .matPun .needGJacExplicit = true; 

or to false respectively. The direct solver for instance sets both values to true, since it 
diagonalizes the Jacobian, the Krylov method sets it so false, since it only needs to evaluate 
the product of the Jacobian with a vector. If the first variable is true, then the global jac 
variable contains the evaluation of the Jacobian at the current step, gjac contains the 
evaluation of gf's Jacobian. 

Additionally, all statistical data fields have to be initialized. They are later used with the 
'statistics' flag. Example: the direct evaluator sets 

eD . stats .matFun. NofDiag = 0; 
eD. stats . matFun. NofMFEv = 0; 

These two variables count the number of matrix function evaluations and the number of 
diagonalizations of the Jacobian. 

We continue with the registerjobs phase. It will be called, after the solver step decides which 
matrix functions it needs to evaluate. This can be in two situations: Either after the matrix 
function evaluator was initialized with the ' init ' flag and now has enough information 
about the differential equation and the options, or when the integration scheme requires a 
change of the matrix function coefficients or even exchanges some matrix functions. The 
latter can be the case for instance in the multistep integrators, after their startup phase is 
finished and the normal integration routine is started. The registerjobs phase is triggered by 

» matFun(jobs, [] , [] , [] , 'registerjobs'); 

The jobs variable is a struct of vectors or matrices. The fieldnames of this struct represent 
the meaning of the vector to multiply the matrix function with. Typically ' F ' is used for a 
product with the right hand side, and ' ^' for the non-autonomous correction. Let job be a 
field of jobs. Each row of the matrix job contains the coefficients for a linear combination of 
the matrix functions in eD . int . jobFunctions. eD . int . jobFunctions contains the scalar, 
vectorized versions of the matrix functions. For exprb it is constructed the following way: 

eD . int . jobFunctions = { 
Ophil 
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@phi2 
@phi3 
@phi4 

@(A, h) phiKA, h/2) 
@(A, h) phi2(A, h/2) 

}; 

The above functions are based on the phim function from the EXP4 package 0. If job is an 
element of jobs for hne number k in a job matrix, matFun should return the following in 
the kth column of the result variable res (written in a suggestive notation): 

length(job{k}) 

res{(:,k)} = ^ job{k}(in) * (eD . int . jobFunctions{m}( J.„ , h) * v) ; 

m=l 

The next flag is the initstep flag: 

>> matFun ([], t, y, h, 'initstep'); 

When this flag is given, the function should prepare itself for several matrix function eva- 
luations with the same t, y and h. 

At the end of the integration process, matFun will be called with flag 'cleanup': 
» matFun ([], [] , [], [] , 'cleanup'); 

Then it should clean up all persistent variables it uses. The matFun field of the global eD 
variable will automatically be cleared by the integrator. 

The call 

>> desc = matFun ([], [] , [] , [] , 'description'); 

should return a basic information string, naming the method used to calculate the results. 
The result will be printed at the beginning of the integration phase by 

>> sprintf ( 'Matrix functions evaluated °/„s.\n', ... 
>> matFun([], [] , [] , [] , 'description')); 

The last flag, which is not used for a computational purpose, is 
»matFun([], [] , [] , [] , 'statistics'); 
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Then matFun is expected to print some statistics. See the ' init flag above for an example 
on what the direct method logs. If there are no interesting statistical data collected, this 
flag should be ignored. 

All other flags that are provided should be used by matFun to recognize the use of the 
executed job. The call looks like this: 

» matFun([], t, y, h, flag, v, reusable, reuse, facs)); 

flag needs to be of the fields of the j obs structure provided in the registerjobs phase. In that 
case, matFun needs to actually compute the result of the matrix functions specified by job 
evaluated at h * J„ or h * A multiplied by the vector v. The required linear combinations 
of the <yC-functions for this computation have already been supplied at the registerjobs phase 
above. 

There are two additional parameters: reusable and reuse. The first one indicates whether 
the matrix function evaluated at the same arguments will have to be used again if the 
current step was rejected by the error estimator. The second one will be set true if the 
previous step was actually rejected and the values calculated there can be reused. In case 
of an exact evaluation, the diagonalization computed the last time can be reused. When 
approximating the results (e.g. in the Krylov version) it may be required to calculate up 
to a higher precision. In case of Kylov approximizations, the old Krylov subspaces can be 
extended. To save the reusable data, use the structure object 

eD .matFun . save . (flag) 

This is one of the reasons to provide those flags. The other reason is to save statistical data 
separately for each matrix function call in the integration scheme. Use 

eD . stats .matFun . (flag) 

as storage here. 

For the addition of exp4 to the EXPODE package, another feature was needed, namely the 
evaluation of ipi{jhJn)v and (p2{jhJn)v, j = 1, facs. This needs to be done at once when 
the facs argument is given with facs > 1. The functions phil and phi2 were extended 
to be called res = phil(v, h, facs) and return an n x 3 Matrix containing the required 
evaluations at once. 

The implementation of custom matrix functions is quite a complex task. The steps needed 
were described as simple as possible even though they are still quite hard to understand 
only by reading this section. Therefore, it is highly recommended to look at the two existing 
implementations as a reference. 
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11 Writing an EXPODE Integrator 



This section contains information for integrator authors. It overviews the EXPODE helper 
routines and how and where to use them. It is split up into several subsections. 

The first subsection describes the fundamentals of the EXPODE package. It names the code 
files needed, explains the package structure and sets coding style standards to make the 
EXPODE code as a whole consistent and well readable. 

The second subsection deals with the integrators themselves. It states the calling conventions 
which will make the usage consistent. It will introduce the helper functions that will ease 
the author's life. These will let him concentrate on the integrator's details instead of having 
to deal with memory management or syntactic correctness of options. Semantic correctness 
still has to be considered by the author. 

The third subsection deals with integrator options. Here we will explain how to create the 
option description structure. 

The last subsection explains the global eD variable. This variable contains all global infor- 
mation shared by EXPODE functions. 

11.1 EXPODE Basics 

There are two different packages for EXPODE. The first one is the normal EXPODE integrator 
package, containing all fully working integrators and all their required helper functions. This 
package is targeted at normal users. The second one, EXPODEDEV is an extension of the first 
one. It additionally contains some tools specifically for developers of new integrators. 

The EXPODE package consists of the following parts: 

• the integrators (expode, exprb, etc), 

• the integrator info functions (exprbinf o, etc), 

• the options helper functions (expset, exprbset, etc), 

• the initPaths routine to setup the appropriate paths, 

• the devalexp function to evaluate the numerical solution at arbitrary times, 

• common helper functions for the integrators in the expode directory, 

• general helper functions in the helpers directory, 

• matrix-function evaluation helpers in the matPun directory, 

• some files from the exp4 package |2] in the Srdparty/ exp4 directory. 
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a set of helper functions for each integrator in folders with the names of the integrators 
and 



• some examples in the examples directory. 

Each integrator has to provide at least the integrator M-file, the options helper function, 
the information command, a function that creates the option description structure and a 
setup command that contains information for the basic integration routine expode and eva- 
luates integrator specific options. The latter two have to be placed in the helper directory 
specific for the integrator. As an example see the files exprb.m, exprbset .m, exprbinf o .m, 
exprb/exprbOpts .m and exprb/exprbSetup .m. The integrator specific helper directory can 
also contain additional supportive functions which will be accessible from the integrator. 

The EXPODEDEV package has the following additional components: 

• a function generateOptionDoku to convert the in-MATLAB help for the options to 
ETgX code for inclusion in a documentation like this one, 

• the BASH script createRelease . sh, which creates release packages, 

• a directory unitTests containing tools for automated tests - useful for finding regres- 
sions, 

• a last directory prototype containing stubs for the minimally required functions for a 
new integrator and a script generatelntegrator . sh to create a new EXPODE integrator 
from them. 

For readability reasons, the code files should obey the following conventions: 

• All M-files should not contain any warnings (and of course no errors) displayed by 
MATLAB's editor. Usually warnings either mark bad coding style or optimization possi- 
bilities. Rarely it is required to disable single warnings like unused function parameters. 

• All names of functions, that are direcly visible to the user should be completely lower 
case. These are mainly the functions in the top level directory except initPaths which 
is used internally to setup paths. Example exprbinf o. 

• All names of variables and internal functions should begin with a lower case let- 
ter. Each new word in the names should start with an upper case letter. Example: 
checkValidOptions, quit Integrator 

• Each M-file should contain a help text in its header. Longer functions that do non- 
trivial work should have comments in the code as well. 

• Indentation is done via four spaces, no tabs. 
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All commas should be followed by a space. Insert spaces after an open squared bracket 
or brace and before a closing one. Empty brackets can be written without spaces. Use 
no spaces around parentheses. Equal signs should be surrounded by spaces. Example: 

» result = str2func([ 'exprb', 'opts', [] ]); 

All lines except the ones that are part of control strucures like if s or whiles should 
be terminated by a semicolon. 

Lines should be maximally around eighty characters log. „Around" means that a line 
should be split up after the first word that exeeds the eightieth column. MATLAB can 
highlight this column. This option can be activated in the preferences dialog under 
„Editor/Debugger" — > „Display". The MATLAB syntax supports split lines via three dots 
at their ends. Strings can be split up, too: 

» text = [ 'This is a very very very ' ... 
» 'very long text' ]; 

Note the space after the last „very" in the first line, because the variable text will not 
contain a newline chacter where the text is split. The remainder of the split lines has 
to be indented eight more spaces than its parent. 

One-line- if constructions should be avoided as well as nested expressions. Generally it 
is dicouraged to put more than one command in a single line. Do not use constructions 
like: 

>> if condition, expression, end 
or 

» varl = funcl(func2(var2, func3(var3) ) , f unc4(var4) ) ; 

The former should be written in three lines, the latter at least expanded into 

» tempi = func2(var2, f unc3(var3) ) ; 

» temp2 = func4(var4); 

» varl = funcKtempl, temp2) ; 

Function code should always be indented one level. 

All functions should be terminated with an end. Nested functions should be indented. 
Additional functions in one M-file should not be intended. The latter ones should only 
be used, if it is clear, that no other function could use the additional function in the 
future. A file majorFunction.m should look like: 
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function majorFunction 
°/o majorFunctionCode 

function nestedFunction 
% nestedFunctionCode 

end 

end 

function additionalFunction 
% additionalFunctionCode 

end 

• There is only one global variable to store common data, a struct called eD, and two 
additional ones to store the Jacobians of the right hand side of the differential, jac 
and a non-linear part in case of seminilinear equations, gjac. EXPODE functions should 
never define additional ones. Generally, it is discouraged to use the eD variable, since 
access to it is slow due to some MATLAB restrictions. Use it only if really necessary or 
the information stored are of common interest doublessly. The two exceptions of this 
only one rule - even though unclean by design - are acceptable since MATLAB has a 
bottleneck accessing big datablocks in global variable structures. 



11.2 Integrators 

In this section we will describe the steps necessary to create a new EXPODE integrator. In 
distinction to the existing integrators we will use a fictitious integrator called expnew here. 

In the EXPODEDEV package there are a number of stub functions for the minimally necessary 
files which need to be created. To create our basic expnew integrator we type 

/path/to/expode $ cd prototype 

/path/to/expode/prototype $ . /generatelntegrator . sh expnew "new-type" 
Its output will look like this: 

Generating integrator expnew, an exponential new-type method. 
[ ... ] 
All done. 

Don't forget to edit the . . /expnew/expnewOpts .m and . ./expnew/expnewSetup.m 
files. Also add the integrator to . ./expode/integratorOpts.m. 

We now have a basic exponential new-type integrator expnew. To activate the integrator, 
edit the file /path/to/expode/expode/integratorOpts .m and add 'expnew' to the list of 
available integrators. Let us first look at the files the script has created. 
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In the /path/to/expode directory we got the new files expnew.m, the integrator command, 
expmsset .m, the options setter and expnewinf o .m, the info command. These files are already 
ready to use, as they only wrap already existing EXPODE functions. These wrappers will 
be adapted by generatelntegrator . sh to give the correct information to their wrapped 
functions. The script gives us a hint where to start our work. All necessary information for 
the integration process have to contained in the two functions in /path/to/expode/expnew, 
expnewSetup .m and expnewOpts .m. The latter one will be discussed in the next section 
where we handle options. 

The expnewSetup. m has to feed the expode in its startup phase where it will be called the 
following way: 

function [ o, nv, solverSetup, solverStep, order, errorOrder, multiStep, . . . 
semilin, denseOutputGenerator ] = expnewSetup (o , nv) ; 

The parameters o and nv contain the options provided by the user and the numeric values 
of the list-type options, see sections 111.31 and 111.41 for details. In case o or nv need to be 
changed by the setup routine, they will be returned again. The arguments solverSetup and 
solverStep will be discussed below in more detail. Next we have order and errorOrder 
which contain the order of the integrator - needed for the user - and the order of the 
error estimator - needed for the step size selection - respectively. Futhermore multiStep 
counts the number of old time steps to use to approximate the solution at the next step, 
semilin indicates, wheter this integrator is built for semilinear problem, expode evaluates 
the linear part once at the beginning instead of evaluating the Jacobian at each time step. 
denseOutputGenerator is a function handle to the dense output generating function. If the 
integrator has no specific one, this can be empty. 

Additionally to these return arguments, the integrator has to set the eD . int . j obFunctions 
field. It has to be a cell array of function handles to the scalar, vectorized versions of the 
matrix functions, that need to be evaluated. The order of the elements has to correspond to 
the order of the job coefficients that are passed to the matrix function evaluator, see section 
[10] for reference. 

It can be necessary to access some of the options set by the user and set some data depending 
on the user's choice. An option's value can be simply accessed via o.OptName. Options of 
tvpe IbooleanI will be evaluated to for false and 1 for true, so you can simply use 

if o . BooleanOption 

°/„ BooleanOption was set 'on' 

else 

y„ BooleanOption was set 'off 

end 

to handle these options. The same works for IlistI type options. You can use the nv variable 
that contains the numeric values corresponding to the list's entries and it is nice to switch 
the option: 
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switch o.List 

case nv. List .value 1 

°/o List was set 'valuel' 
case nv.List . value2 

% List was set value2' 
case nv.List .values 

°/„ List was set 'value 3' 

end 

We now discuss the missing two return arguments solverSetup and solverStep of the 
setup command, which both are function handles. The first one will be called to initialize 
the solver step once, after all initialization is complete and the integration process is about 
to begin, just before the initialization of the matrix function evaluator. This function gets 
no arguments. For simple, static integrators, where there is no need for such a function, it 
can be set to Snullf unction. The solverStep performs the actual integration step. It has 
to have the following signature: 

function [ yNew, normErr, hOut ] = solverStep (t , y, h, reuse) 

Here we obviously have the time t and phase variable y, where we start our step with step 
width h. The output yNew is the numerical solution at time t + hOut. Note that it is possible 
that the step size had to be reduced from h to hOut during the computation of the new step, 
for instance if a Krylov method was used. The output argument normErr contains the norm 
computed by the error estimator, if one is in use, set to otherwise. The reuse argument 
is true if the preceeding step was rejected due to the error estimator. It should be passed 
to the matrix function evaluator for all complutations that are only based on the old time 
step. In case it is true some of the previously computed could be reused. 

It is generally a good idea to have a look at the existing implementations of solver steps, 
such as exprb32, ms or exp4, to get a good impression how they have to work. The use of 
the matrix function evaluators was already explained in section [TOl Note that it might be 
possible that the matrix function evaluator requires to reduce the step size. In this case all 
previously computed data for the next time step that is influenced by the step size have to 
be recalculated. That is the reason for the while-loops in all the existing solver steps. 

You can access old and the current evaluations of the right hand side and possibly of the 
non-linear part via 

» [ Fl, Gl ] = getOldF(O); 

where the argument of getOldF determines how many steps you need to go to the past, 
starting with for -F(?/„) only. 

Now you should have the basic information you need to write expnewStep. The remainder 
of this document deals with the definition of user options and the usage of the global eD 
variable which are the two crutial parts left to understand to complete the integrator. 
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11.3 Options 



This section deals with integrator options. EXPODE has advanced option handhng subroutines 
which allow automatic checking of options set by the user. First, each integrator needs an 
options helper function as a user interface. This function should be called the same as the 
integrator with a „set" suffix. The exprb integrator is contained in the exprb . m M-file and 
its options helper is named „exprbset .m" for instance. The options helper should simply 
be a wrapper around the expset function which only adds the integrator name as the first 
argument : 

function opt = expnewset (varargin) 
intname = 'expnew'; 
initPaths (intname) ; 

if isempty (varargin) 
help expnewset ; 
return; 

end 

opt = expset (intname , varargin{ : }) ; 

end 

The initPaths routine is needed to initialize the EXPODE internal paths as seen in the pre- 
vious section. This method can be automatically generated by the generatelntegrator . sh 
script, see the previous section. 

The expset method needs a so-called option description object. This has to be provided by 
a function called expnewOpts and has to be located in the expnew directory. This object 
contains all information about the options for expnew. It is a MATLAB struct object and has 
to have the following strutcure: 

>> expnewOpts 



The first four fields are all strings, the first one containing the integrator name, the second 
a longer name of the integration method, the third a short description, and the last one 
the help text for the integrator, displayed when calling the integrator without arguments. 



ans = 



integrator : 
integratorname : 
integratordesc : 
usage : 
opts : 



' expnew ' 

'exponential new-type integrator' 



[1x223 char] 
[1x522 char] 
[1x1 struct] 
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This will be done automatically when wrapping the expode base integrator, see the previous 
section. 

The last entry, opts, is again a structure object and it contains the actual options. Each 
element of this structure represents an option where the name of the element ist the name 
of the option. The option itself is a struct once again and can be created via 



opts . NewOpt ion = createOptDesc( ... 
'short descripuion , ... 
'type', ... 
'default', ... 

{ 'valuel' , 'value2' }, { ... 

'Very long and detailed description for the very nice' 
'and new option NewOpt ion.' 
}. { 

'OptionA', 'OptionB' 



The first argument is a short description that will be displayed when using expnewinf o, see 
below. The next one contains the option's type. Available option types are listed in section 
Ei Next, we have the default value for the option. This can be a numeric value or string, 
depending on the option type. In case of Ibooleaiil or IlistI it can either be the numeric 
or string value of the option. Argument four of the createOptDesc function contains the 
allowed values for IlistI For all other types, provide the empty cell object {}. The next 
argument contains the long and detailed description for the option. This can either be a 
string or a cell of strings as used in the example. The last argument contains a list of 
referenced options. This list will be printed in the form See also: OptionA, OptionB 
when calling expnewinf o NewOpt ion. 

For advances usage two additional arguments can be given, renameTo and aliases, in that 
order. The first one is used at the moment, where options are evaluated to their numeric 
values. The option will then be renamed to renameTo. This is useful if the same option has 
different names in different scenarios. In case of variable step size integrators, there is an 



option [Tnit ialStepI which indicates the step size to start with, where in the case of constant 



step size options there is the |StepSize| option, which will be renamed to InitialStep in 
the integration process to avoid code duplication. 

For reference and a lot of examples look at the file expode/baseOpts .m or any other 
*/*Opts.m file in the EXPODE package. 

A basic stub for the opts command will be generated by the generatelntegrator . sh script. 
There are a set of basic options the integrator has to provide, see section IHfor details. They 
can be generated with 



optionDesc = baseOpts( 'expnew' , 'exponential new-type integrator'); 
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This line of expnewOpts .m is part of the stub and doesn't need to be adjusted. These basic 
options have to be extended either for use in semihnear integrators or hnearized integrators 
and for constant or variable step size integrators. The two code lines responsible for this are 

optionDesc = linearizationOpts (optionDesc) ; 
optionDesc = varStepOpts (optionDesc ) ; 

Use semilinOpts instead of linearizationOpts for semilinear integrators, constStepOpts 
instead of varStepOpts for constant step size integration. Using constStepOpts automati- 
cally switches off the step size estimator. 

Another function that needs to be provided for each integrator is the info command. This 
is a wrapper around optinf o similar to the options setter. It will print information on the 
integrator. The information will be extracted from the option description object described 
above. Here is its basic structure: 

function expnewinf o(optName) 
if nargin < 1 

optname = ; 

end 

optinf o( 'expnew' , optName) ; 

end 

Like the set, opts and the integrator command this file will be generated by the script 
generatelntegrator . sh. 

11.4 The Global eD Variable 

This section deals with the main global variable eD. Except for jac and jacv this is the 
only global variable used throughout EXPODE (see the coding conventions above). Since it 
has such basic relevance it will be discussed in detail here. It is a structure object containing 
several different types of information. After an integrator call, here exprb as an example, it 
contains the following fields if the IClearlnternalDatal is set to ' of f ' : 



» eD 
eD = 



integrator 


' exprb ' 


ODE 


[1x1 


struct] 


int 


[1x1 


struct] 


log 


[1x1 


struct] 


evals 


[1x1 


struct] 


functions 


[1x1 


struct] 
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stats: [1x1 struct] 

matFun: [1x1 struct] 

tdata: [451x1 double] 

ydata: [451x19 double] 

The first field is integrator and it holds the name of the integrator. The second, ODE, 
contains information about the differential equation and the solver request: 

» eD.ODE 



ans = 



to 





T 


45 


duration 


45 


tDir 


1 


yo 


[19x1 double] 


ylen 


19 


transpose 


1 


options 


[1x1 struct] 



It contains the endpoints of the integration interval, tO and T, their difference duration, the 
integration direction (tDir) - with values 1 for increasing or -1 for decreasing time steps, the 
initial value yO and its length ylen as well as the options provided at the integrator's call, 
transpose is a boolean that indicates, whether the solution vectors are rows or columns. In 
the latter case, they have to be transposed before being stored in the result matrix. 

The int field of eD contains information for the integration process. Is has the following 
form: 

» eD.int 

o : 
nv: 
od: 

jobFunctions : 
order : 
errorOrder : 
mult i Step : 
semilin : 
denseOutputGenerator : 
wantSolverOutput : 
wantDenseOutput : 
solverOutputStepData : 



[1x1 struct] 

[1x1 struct] 

[1x1 struct] 

6x1 cell 

4 

4 

1 



Ohermite 



1 

{[19x1 double] [19x1 double]} 
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The fields o contains a version of tlie options structure but witli all non-set options evaluated 
to their default values and all llistl -tvpe options evaluated to the index of the selected value 
in the list. Ibooleaiil -type options are also converted to 1 for value true and for false. The 
last field, nv, contains the numeric values of the list-type options, see the previous section. 
The field od contains the option description strcture describing the integrator's options 
that were used to evaluate eD. ODE. options into eD.int.o. jobFunctions contains function 
handles to scalar, vectorized versions of the matrix functions involved in the integrator's 
scheme. These will be used by matrix function evaluators, see sections [10] and and II 1.21 The 
field order represents the integrator order, whereas errorOrder ist the order of the error 
estimator plus one - i.e. the parameter for the errstep function used to compute the next 
timestep size. The entry multiStep contains the number of old time step data necessary to 
compute the next time step, semilin indicates, whether the solver is for semilinear problems, 
i.e. whether in needs to evaluate the linear part once at the beginning or the Jacobian at 
each time step. The denseOutputGenerator field is a handle to the dense output generation 
function. The two flags wantSolverOutput and wantDenseOutput indicate, wheter the user 
requested one of the two outout types. If the first is true, then solverOutputStepData 
contains the data - additionally to the time step data - needed to evaluate the numerical 
solution by the devalexp function. If the second is true it contains the data necessary for 
the dense output evaluation between the last two time steps. 

Next in the eD structure we have the log field containing logging functions for several 
purposes. They all have the same syntax as the f printf command, but may be directed to 
a file or simply ignored depending on what the user wants to see. 



» eD.log 
ans = 

verbose : 
status : 
statistics : 
jacLog: 
stepLog : 
matFunLog : 
warning : 
error : 



Of printf 

(varargin) wrapped (fun,prefix,varargin{ : }) 

Of printf 

Onullf unction 

Of printf 

Onullf unction 

Owarning 

Oerror 



In eD. functions we can find evaluation functions for the ODE and its Jacobian, furthormore 
we find a handle to the matrix function evauator and last we can retrieve the varargin 
parameters provided by the user at the integrator's call. In the example the field reads 



>> eD. functions 
ans = 

ode : (t , y , varargin) ode (t , y , flag , varargin-[ : }) 
jac : 0(t ,y, varargin) ode (t ,y, flag, varargin-[ : }) 
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matFunV : SmatFunKrylov 
varargin : {} 



The stats field in eD contains statistical data collected over the whole integration process. 
It has a subfield matPun which contains data gathered by the matrix function evaluators. 

The fields tdata and ydata contain the data of the numerical solution. This is what the 
integrator returns to the user. 

Most of this data will be cleaned before the integration ends. If you want to preserve the 
data for testing, set the option IClearlnternalDatal to . 

11.5 Deploying EXPODE 

The previous subsection concludes the basic information needed to write new integrators. 
Additional help can be found using MATLAB's help command on the EXPODE functions. All 
should be documented well both helptext and code wise. In this subsection you will find 
some advices for deployment of your newly created integrator. 

The EXPODEDEV package contains a unit testing framework. Before packaging the new release 
it is highly recommended that you run at least the existing tests to make sure that you didn't 
introduce any regression to the previous release. This is done by 

» cd unitTests; 
» runAllTests; 

If all goes well you should see a message like 
Expode unit tests framework finished. 

Totally ran 593 tests, out of where 593 finished successful and failed. 
Overall testing time: 00:01:37. 

Additionally you should write new unit tests for you new components und run them. Tests 
have to be located in the unitTests/tests directory. The testing framework clears all 
global variables and restores MATLAB's default path before calling the actual test. Only one 
global variable for the testing data will be created, the test itself does not have to use it, 
though, there are helper functions for it. The test's working directory is unitTests/tests. 
You will need to call an intialization routine to tell the framework what you want to test. 
The next thing you need to do is add the path to your tested functions. 

startTest ( ' optionDesc ' ) ; 
addpath . . / . . /expode ; 

Now you can write your testing code. You may need set some values in the global eD variable. 
After you are done, use 
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finishTest (success) ; 



to report the result, where success is a boolean containing your result. Running this specific 
test is done via 

>> runTest testname; 

run from the unitTests directory again, where testname. m is the name of the m-file. Note 
that you can include more that one test into one file, you just have to separate them by 
startTest and finishTest. The runTest routine can also run multiple tests in a row, just 
give more testnames separated by a space. 

Now that you finished testing you new integrator, it is time to release it. In EXPODEDEV you 
will find a script, createRelease . sh to do that. You have to add the integrator to the ints 
variable. This will automatically include the integrator m-file, the options helper, the info 
command and the integrator's directory. If you have additional files, you have to add it to 
the apropiate variable in the script. This should usually not be the case, since all methods 
that are only used by your new integrator should be in its own directory. Now you can call 

/path/to/expode $ . /createRelease . sh 1.1 

to create a package for EXPODE, version 1.1. To generate a EXPODEDEV package, call 
/path/to/expode $ . /createRelease . sh 1.1-dev -d 

which will result in version 1.1-dev with the additional functions described at the beginning 
of this section. 
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